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Developing a theory of low-mass star formation (~ 0.1 to 3 M0) remains one of the 
' most elusive and important goals of theoretical astrophysics. The star-formation process is 

the outcome of the complex dynamics of interstellar gas involving non-linear interactions of 
turbulence, gravity, magnetic field and radiation. The evolution of protostellar condensations, 
from the moment they are assembled by turbulent flows to the time they reach stellar densities, 
spans an enormous range of scales, resulting in a major computational challenge for simulations. 
Since the previous Protostars and Planets conference, dramatic advances in the development 
, of new numerical algorithmic techniques have been successfully implemented on large scale 

■ parallel supercomputers. Among such techniques. Adaptive Mesh Refinement and Smooth 
' Particle Hydrodynamics have provided frameworks to simulate the process of low-mass star 
, formation with a very large dynamic range. It is now feasible to explore the turbulent fragmen- 
' tation of molecular clouds and the gravitational collapse of cores into stars self-consistently 

within the same calculation. The increased sophistication of these powerful methods comes 
^ ^. with substantial caveats associated with the use of the techniques and the interpretation of the 

I ' numerical results. In this review, we examine what has been accomplished in the field and 

O I present a critique of both numerical methods and scientific results. We stress that computational 

simulations should obey the available observational constraints and demonstrate numerical 



^ ' convergence. Failing this, results of numerical simulations do not advance our understanding of 

. . , low-mass star formation. 

> . 

■ 1. INTRODUCTION inherent in the star formation process have made it difficult 

to perform accurate calculations of fragmentation and col- 
Most of the stars m the galaxy exist in eravitationally , , . , ■ . • • n j- • 1 • 

° -' ° -' lapse, which are intrinsically three-dimensional in nature, 

bound binary and low-order multiple systems. Although „. 1 . • • n j m tt/i. r. 

■' f J o Since the last review in Protostars and Planets IV by Bo- 
several mechanisms have been put forth to account for bi- , , • , , ^nnn j j • .i, j 1 

^ denheimer et ai, 2000, dramatic advances in the deve lop- 
nary star formation, fragmentation has emerged as the lead- ^ J. 11 -.u • . u • ■ 1 J- 

■' ' o o ment of new numerical algorithmic techniques, including 

mg mechanism in the past decade (Bodenheimer et ai, , . . ^ waajon j c .u n 1 tt 

±- V adaptive mesh rehnement (AMR) and Smooth Particle Hy- 

2000). This point of view has been strengthened by obser- . . . /cnii\ n j 1 1 

^ b J drodynamics (SPH), as well as advances in large scale par- 

vations that have shown the binary frequency among pre- „ , , . , n j • vc . • • 

^ , , ^llcl machines, have allowed a signmcant increase in the 



main-sequence stars is comparable to or greater than that 
among nearby main-sequence stars (Duchene et ai, 1999). 
This suggests that most binary stars are formed during the 
protostellar collapse phase. Developing a theory for low 



dynamic range of simulations of low mass-star formation. 

It is now feasible to explore the turbulent fragmentation of 

molecular clouds and the gravitational collapse of cores into 

stars self-consistently within the same calculation. In this 
mass star formation ( 0. 1 to 3 solar masses), which explains . u.uu »i i-uj 

^ r chapter we examine what has been recently accomplished 
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. t- J f f J jjjg of numerical simulation of low-mass star forma- 

tiple stellar systems, remains one of the most elusive and . , n • 1. • 1 .u j j 

^ ■' tion, and we critically review both numerical methods and 

important goals of theoretical astrophysics. . ,^ 

scientitic results 

Until very recently, the extreme variations in length scale 
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1.1. Key Questions Posed by the Observations 

Observational surveys present us with a basic picture 
of star-forming regions including the structure and dynam- 
ics of star-forming clouds and the properties of protostellar 
cores. A theory of star-formation should explain both the 
large scale environment and the properties of protostellar 
cores self-consistently. 

As shown first by Larson (1981), and later confirmed by 
many other studies, star-forming regions are characterized 
by a correlation between internal velocity dispersion and 
size, 5V « lkm/s(L/lpc)°-^. This scaling law has been 
interpreted as evidence of supersonic turbulence on a wide 
range of scales. The turbulence can provide support against 
the gravitational collapse, but can also create gravitationally 
unstable compressed regions through shocks. A theory of 
star formation should elucidate whether turbulence controls 
the star formation rate and efficiency, or are those properties 
controlled by stellar outflows and winds? Are the scaling 
laws of turbulent flows related to scaling laws of core and 
stellar properties? Is the turbulence setting the initial den- 
sity perturbations that collapse into stars? What is the effect 
of turbulence on the accretion of mass onto protostars? 

On smaller scales, observational surveys have shown 
that prestellar cores are elongated. Their density profiles are 
flat near the center, steeper at larger radii, may show very 
sharp edges and are sometimes consistent with Bonnor- 
Ebert profiles (e.g., Bacmann et al, 2000; Alves et al., 
2001). Cores have rotational energies on the average only 
a few percent of their binding energies (e.g., Goodman et 
al, 1993). They are marginally supercritical {Cmtcher, 
1999) and their mass distribution is very similar to the stel- 
lar mass distribution (e.g., Motte et al, 1998, 2001). The 
large majority of cores are found to contain stars (e.g., Ji- 
jina et al, 1999), and individual cores produce at most 2-3 
protostars (e.g., Goodwin and Kroupa, 2005). A theory of 
low-mass star formation should be consistent with these ob- 
servations and address: Why are prestellar cores so short- 
hved? Why do they have Bonnor-Ebert profiles? How does 
the observed core angular momentum affect the formation 
of binary and multiple systems? What is the role played by 
magnetic fields in their evolution? Why is the mass distribu- 
tion of prestellar cores so similar to the stellar initial mass 
function? Why are cores barely fragmenting into binaries 
or low multiphcity systems? 

Observations of young stellar populations provide im- 
portant constraints as well. We know that young stars are 
always found in association with dense gas with an effi- 
ciency ~ 10-20%, much higher than the overall star forma- 
tion efficiency in GMCs, ^ 1-3% (e.g., Myers et al, 1986). 
Stars are often found in clusters, of size ranging from 10 
to 1000 members (e.g., Lada and Lada, 2003). The stel- 
lar initial mass function peaks around a fraction of a solar 
mass and its lognormal shape around the peak is roughly 
the same in open clusters, globular clusters and field stars 
(Chabrier, 2003). What determines the efficiency of star 
formation? Why is the stellar multiplicity higher in younger 



populations? What determines the typical stellar mass and 
the initial mass function? 

Although answering all these questions is outside the 
scope of this review, we pose them because these are the 
questions to be addressed by the computational simulations 
of the formation of low-mass stars. 

1.2. Generation of Initial Conditions Consistent with the 
Observations 

Simulations of low mass star formation should gener- 
ate initial conditions for the collapse of protostars consis- 
tent with the observed physical properties of star forming 
clouds. This can be achieved if a relatively large scale is 
simulated (> 1 pc) with numerical methods that can accu- 
rately reproduce fundamental statistics measured in molec- 
ular clouds. Such statistics include i) scaling laws of veloc- 
ity, density and magnetic fields; ii) mean relative values of 
turbulent, thermal, magnetic and gravitational energies (the 
normalization of the scaling laws). 

7.2.7. Scaling Laws. Larson (1981) found that veloc- 
ity and size of interstellar clouds are correlated over many 
orders of magnitude in size. This correlation has been con- 
firmed by many more recent studies (e.g.. Fuller and Myers, 
1992; Falgarone et al, 1992). The most accepted interpre- 
tation is that the scaling law reflects the presence of super- 
sonic turbulence in the ISM (e.g., Larson, 1981; Ossenkopf 
and Mac-Low, 2002; Heyer and Brunt, 2004). Large scale 
velocity-column density correlations from molecular line 
surveys of giant molecular clouds also suggest a turbulent 
origin of the observed density enhancements (Padoan et al , 
2001). Starting with the work of Troland and He iles (1986), 
a correlation between magnetic field strength and gas den- 
sity, B oc n^/^, has been reported for mean densities larger 
than n 100 cm~'^. However, density and magnetic field 
scalings are very uncertain because both quantities are dif- 
ficult to measure. 

7.2.2. Mean Energies. Assuming an average tempera- 
ture of T = 10 K, Larson's velocity-size correlation cor- 
responds to an rms sonic Mach number Mg « 5 on the 
scale of 1 pc, and Mg « 1 at 0.02 pc. So, on the average, 
the turbulent kinetic energy is larger than the thermal en- 
ergy. Indirect evidence of super- Alfvenic dynamics in giant 
molecular clouds has been presented by Padoan and Nord- 
lund (1999) and Padoan et al. (2004a). They have shown 
that the magnetic energy averaged over a large scale has an 
intermediate value between the thermal and the kinetic en- 
ergies, even if it can be significantly larger than this average 
value within dense prestellar cores. Observations suggest 
that in dense cores gravitational, kinetic, thermal and mag- 
netic energies are all comparable. However, the magnetic 
energy is very difficult to estimate. Accounting for both 
detections and upper limits, there is a large dispersion in 
the ratio of magnetic to gravitational energy of dense cores 
(Crutcher et al, 1993; Crutcher et al, 1999; Bourke et al, 
2001; Nutter et al, 2004). In the case of turbulence that 
is super- Alfvenic on the large scale, this dispersion and the 
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B-n relation are predicted to be real {Padoan andNordlund, 
1999). 

In summary, large scale simulations may provide real- 
istic initial and boundary conditions for protostellar col- 
lapse, but they must be consistent with the turbulent nature 
of the ISM. On the scale of giant molecular clouds, obser- 
vations suggest the turbulence is on the average supersonic 
and super-Alfvenic and its kinetic energy is roughly equal 
to the cloud gravitational energy. 

2. A BRIEF SURVEY OF LOW-MASS STAR FOR- 
MATION MECHANISMS 

Although much progress in numerical simulations of 
collapse and fragmentation has been made in the interven- 
ing 6 years since PPIV, a self-consistent theory of binary 
and multiple star formation that addresses the key questions 
posed by observations is still not at hand. As discussed 
by Bodenheimer et at, 2000, binary and multiple forma- 
tion can occur through the processes of (i) capture, (ii) fis- 
sion, (iii) prompt initial fragmentation, (iv) disk fragmenta- 
tion and (v) fragmentation during the protostellar collapse 
phase. 

A recent mechanism for multiple star formation has been 
put forth by Shu et al, (2000) and Galli et al, (2001). 
They develop equilibrium models of strongly magnetized 
isopedic disks and explored their bifurcation to non-axi- 
symmetric, multi-lobed structures of increasing rotation 
rates. Possible problems with this mechanism include 
the observed low rotational energies, the observed ran- 
dom alignment of disks with the ambient magnetic fields, 
the complexity of star-forming regions relative to the two 
dimensional geometry and absence of turbulence in the 
model. 

Disk fragmentation from gravitational instability can re- 
sult in multiple systems in an equiUbrium disk if the mini- 
mum Toomre Q parameter falls below « 1. However, Bo- 
denheimer et al, 2000, have pointed out that the required 
initial conditions to obtain Q <1 may not be easily reaUzed 
since the mass accretion timescale is significantly longer 
than the dynamical timescale throughout most of the evo- 
lution of the protostar. Disk fragmentation plays a key role 
in one of the theories of the formation of Brown Dwarfs 
(BDs). This scenario, known as the "failed embryo" sce- 
nario, begins with a gravitationally unstable disk surround- 
ing a protostar. The disk fragments into a number of sub- 
stellar objects. If the crossing time of the cluster of embryos 
is much less than the free-fall time of the collapsing core, 
one or more of the members will be rapidly ejected result- 
ing in a BD (Reipurth and Clarke, 2001). Problems with 
this model include observational evidence of BD cluster- 
ing {Duchene et al., 2004) , Ly-a signatures of BD accre- 
tion (e.g., Jayawardhana et al., 2002; Natta et al., 2004; 
Barrado y Navascues et al, 2004; Mohanty et al., 2005) 
and evidence that individual cores produce only 2 or 3 stars 
{Goodwin and Kroupa, 2005). 

Currently, there are two dominant models to explain 



what determines the mass of stars. The Direct Gravitational 
Collapse theory suggests that star-forming turbulent clouds 
fragment into cores that eventually collapse to make indi- 
vidual stars or small multiple systems {Shu et al, 1987; 
Padoan and Nordlund, 2002, 2004). In contrast, the Com- 
petitive Accretion theory suggests that at birth all stars are 
much smaller than the typical stellar mass and that the fi- 
nal stellar mass is determined by the subsequent accretion 
of unbound gas from the clump {Bonnell et a/., 1998; Bon- 
nell et al., 2001; Bate et al, 2005). Significant problems 
with competitive accretion include the large value of the ob- 
served virial parameter relative to the one required by com- 
petitive accretion {Krumholz etal, 2005b). We discuss this 
problem with competitive accretion in detail in section 4d. 

3. PHYSICAL PROCESSES NECESSARY FOR 
DETAILED SIMULATION OF LOW-MASS STAR 
FORMATION 

3.1. Turbulence 

The Reynolds number estimates the relative importance 
of the nonlinear advection term and the viscosity term in 
the Navier-Stokes equation. Re = VqLq/v. Vq is the flow 
rms velocity, Lq is the integral scale of the turbulence (say 
the energy injection scale) and v is the kinematic viscos- 
ity that we can approximate as « Vth/{<7n). Vth is the 
gas thermal velocity, n is the gas mean number density and 
(T ^ 10^^^ cm^ is the typical gas collisional cross section. 
For typical molecular cloud values. Re ~ 10''-10®, which 
implies flows are highly unstable to turbulence. It is impor- 
tant to recognize the significance of turbulent gas dynam- 
ics in astrophysical processes, as turbulence is a dominant 
transport mechanism. In molecular clouds, the turbulence is 
supersonic and the postshock gas cooling time is very short. 
This results in the highly fragmented structure of molecular 
clouds. 

There has been significant progress in our understand- 
ing of supersonic turbulence in recent years (progress on 
the scaling properties of subsonic and sub-Alfvenic turbu- 
lence is reviewed in Cho et al, 2003). Phenomenological 
models of the intermittency of incompressible turbulence 
(e.g.. She and Leveque, 1994; Dubrulle, 1994) have been 
extended to supersonic turbulence by Boldyrev (2002) and 
the predictions of the model have been confirmed by nu- 
merical simulations {Boldyrev et al, 2002a; Padoan et al, 
2004b). The intermittency correction is small for the expo- 
nent of the velocity power spectrum (corresponding to the 
second order velocity structure function) and large only at 
high order. However, Boldyrev et al. (2002b) have shown 
that low order density correlators depend on high order ve- 
locity statistics, so intermittency is likely to play a signif- 
icant role in turbulent fragmentation, despite being only a 
small effect in the velocity power spectrum. 

Because supersonic turbulence can naturally generate, 
at very small scale, strong density enhancements of mass 
comparable to a low mass stars or even a brown dwarf. 
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its correct description is of paramount importance for sim- 
ulations of molecular cloud fragmentation into low mass 
stars and brown dwarfs. At present, the largest simula- 
tions of supersonic turbulence may achieve a Reynolds 
number Re ~ 10"^. The scale of turbulence dissipa- 
tion is therefore much larger in numerical simulations 
(of order the computational mesh size) than in nature 
(~ 10^^ cm). However, the ratio of the Kolmogorov dis- 
sipation scale,77K, and the Jeans length, Aj, is very small 
and remarkably independent of temperature and density, 
tjk/Xj « 10-''(r/10K)-i/8(n/103cm-3)-i/4, One may 

hope to successfully simulate the process of turbulent frag- 
mentation by numerically resolving the turbulence to scales 
smaller than Aj, but not as small as rjK, unless the nature of 
turbulent flows varies dramatically between Re ^ 10^ and 
Re ~ 10^. Experimental results seem to indicate that the 
asymptotic behavior of turbulence is recovered in the ap- 
proximate range Re = 10, 000-20,000 (Dimotakis, 2000), 
which can be achieved with PPM simulations on a 2, 048^ 
mesh. 

In order to i) generate a sizable inertial range (a power 
law power spectrum of the turbulent velocity over an ex- 
tended range of scales) and ii) resolve the turbulence just 
below the Jeans length, a minimum computational box size 
of at least 1, 000^ zones is required for a grid code. This ac- 
counts for the fact that the velocity power spectrum starts to 
decay with increasing wavenumber faster than a power law 
already at approximately 30 times the Nyquist frequency. 
It is still an optimistic estimate, because at this resolution 
the velocity power spectrum is also affected by the bot- 
tleneck effect (e.g., Falkovich, 1994; Dobler et al, 2003; 
Haugen and Brandenburg, 2004). Assuming the standard 
SPH kernel of 50 particles, this corresponds to at least 
50 X 1, 000^ particles to describe the density field, and at 
least few x 1 , 000^ particles to describe the velocity field, 
if a Godunov SPH method is used (see below). 

Grid code simulations have started to achieve this dy- 
namical range only recently, while particle codes appear un- 
suitable to the task. The calculation of Bate et al. (2003) 
has 3.5 X 100^ particles, more than four orders of magni- 
tude below the above estimate and therefore inadequate to 
describe the process of turbulent fragmentation (the forma- 
tion of small scale density enhancements by the supersonic 
turbulence). Studies proposing to directly test the effect of 
turbulence on star formation, based on numerical simula- 
tions with resolution well below the above estimate, should 
be regarded with caution. 

3.2. Gravity 

3.2.1. The Jeans Condition. Jeans (1902) analyzed the lin- 
earized equations of ID isothermal self-gravitational hydro- 
dynamics (GHD) for a medium of infinite extent and found 
that perturbations on scales larger than the Jeans length, 
Aj = (^^^)^^^> are unstable. Thermal pressure cannot 
resist the self-gravity of a perturbation larger than Aj, re- 
sulting in runaway collapse. Truelove et al. (1997) showed 



that the errors generated by numerical GHD solvers can act 
as unstable perturbations to the flow. In a simulation with 
variable resolution, ceU-scale errors introduced in regions 
of coarser resolution can be advected to regions of finer res- 
olution, allowing these errors to grow. The unstable col- 
lapse of numerical perturbations can lead to artificial frag- 
mentation. The strategy to avoid artificial fragmentation is 
to maintain a sufficient resolution of Aj. Defining the Jeans 
number J = Truelove et al. (1997) found that keeping 
J < 0.25 avoided artificial fragmentation in the isothermal 
evolution of a collapse spanning 7 decades of density, the 
approximate range separating typical molecular cloud cores 
from nonisothermal protostellar fragments. This Jeans con- 
dition arises because perturbations on scales above Aj are 
physically unstable, and discretization of the GHD PDEs 
introduces perturbations on all scales above Ax. It is es- 
sential to keep Aj as resolved as possible in order to dimin- 
ish the initial amphtude of perturbations that exceed this 
scale. Although it has been shown to hold only for isother- 
mal evolution, it is reasonable to expect that it is necessary 
(although not necessarily sufficient) for nonisothermal col- 
lapse as well where the transition to non-isothermal evolu- 
tion may produce structure on smaller scales than the local 
Jeans length. . 

3.2.2. Runaway Collapse. The self-gravitational collapse in 
nearly spherical geometry tends to show a so-called "run- 
away collapse,", where the denser central region collapses 
much faster than the less-dense surrounding region. The 
mass of the central fast collapsing region is of the order of 
the Jeans mass, Mj = pA^ ~ G-^/^cV^ p-^'^, which de- 
creases monotonically in this runaway stage. The descrip- 
tion of this process requires increasingly higher resolution, 
not only on the spatial scale but also on the mass scale. 
Therefore, an accurate description is not guaranteed even 
with Lagrangian particle methods such as SPH, if the num- 
ber of particles is conserved. The end of the runaway stage 
corresponds to the deceleration of the gravitational collapse. 
If the effective ratio of specific heats, 7 (P ex p^), becomes 
larger than 7crit = 4/3, the increased pressure can decel- 
erate the gravitational collapse. For example, the question 
of how and when the isothermal evolution terminates was 
explored in Masunaga andlnutsuka (1999). 

3.2.3. Thermal Budget. In the low density regime the gas 
temperature is affected by various heating and cooling pro- 
cesses (e.g., Wolfire et al., 1995; Koyama and Inutsuka, 
2000; Juvela et al., 2001). However, above a gas density 
of 10^-10^ cm~"^, depending on the timescale of interest, 
the gas is thermally well coupled with the dust grains that 
maintain a temperature of order lOK. During the dynami- 
cal collapse, gas and dust are isothermal until a density of 
10^°-10^^ cm~'^, when the compressional heating rate be- 
comes larger than the cooling rate (Inutsuka and Miyama, 
1997; Masunaga and Inutsuka, 1999). The further evolu- 
tion of a collapsing core and the formation of a protostar are 
radiation-hydrodynamical (RHD) processes that should be 
modeled by solving the radiation transfer and the hydrody- 
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namics simultaneously and in three dimensions. F*resently, 
the most sophisticated multi-dimensional models are based 
on the (flux-hmited) diffusion approximation (Bodenheimer 
etal. 1990, Kmmholz et al. 2005c). 

Once the compressional heating dominates the radiative 
cooling, the central temperature increases gradually from 
the initial value of ~ 10 K. The initial slope of the tem- 
perature as a function of gas density corresponds to an ef- 
fective ratio of specific heats 7 = 5/3: T{p) cx p^/^ for 
lOK < T < lOOK. This monatomic gas property is due 
to the fact that the rotational degree of freedom of molecu- 
lar hydrogen is not excited in this low temperature regime 
(e.g., E{J = 2 - 0)/fcB = 512K). When the tempera- 
ture becomes larger than ^ lO^K, the slope corresponds 
to 7 = 7/5, as for diatomic molecules. This value of 7 is 
larger than the minimum required for thermal pressure sup- 
port against gravitational collapse: 7 > 7crit = 4/3. The 
collapse is therefore decelerated and a shock is formed at 
the surface of a quasi-adiabatic core, called "the first core". 
Its radius is about 1 AU in spherically symmetric mod- 
els, but can be significantly larger in more reahstic multi- 
dimensional models. It consists mainly of H2. 

The increase of density and temperature inside the first 
core is slow but monotonic. When the temperature becomes 
> 10^ K, the dissociation of H2 starts. The dissociation 
of H2 acts as an efficient cooling of the gas, which makes 
7 < 4/3, triggering the second dynamical collapse. In 
this second collapse phase, the collapsing velocity becomes 
very large and engulfs the first core. As a result, the first 
core lasts for only ~ 10'^ years. In the course of the second 
collapse, the central density attains the stellar value, p ~ 
1 g/cm^, and the second adiabatic core, or "protostar", is 
formed. The time evolution of the SED obtained from the 
self-consistent RHD calculation can be found in Masunaga 
et al. (1998) and Masunaga and Inutsuka (2000a,b). 
3.2.4. Core Fragmentation. Tsuribe and Inutsuka (I999a,b) 
have shown that the fragmentation of a rotating collaps- 
ing core into a multiple system is difficult in the isother- 
mal stage. Matsumoto and Hanawa (2003) have extended 
the collapse calculation by using a nested-grid hydro code 
and a simplified barotropic equation of state that mimics the 
thermal evolution, and have shown that the first-core disk 
increases the rotation-to-gravitational energy ratio (T/ \W\) 
by mass accretion. A stability analysis of a rotating poly- 
tropic gas shows that gas with T/|M^| > 0.27 is unstable 
for non-axisymmetric perturbations (e.g., Imamura et al. 
2000). If the first-core disk rotates fast enough that the an- 
gular speed X the free-fall time satisfies ric(47rGpc)~^^^ ^ 
(0.2 — 0.3), fragments appear and grow into binaries and 
multiples in the first core phase. The non-axisymmetric 
nonhnear spiral pattern can transfer the angular momentum 
of the accreting gas. 

3.3. Magnetic Fields 

Detailed self-consistent calculations accounting for 
both thermal and magnetic support (Mouschovias and 



Spitzer, 1976; Tomisaka et al, 1988) show that the max- 
imum stable mass can be expressed as Mmag.max ~ 

{9 —3/2 
1 - [0.17/(GV2m/$)c] I , where (M/$)c is 

the central mass-to-flux ratio and Mbe = l-lSCg/G'^^^Pg^t 
is the Bonnor-Ebert mass (Bonnor, 1956, 1957; Ebert, 
1957). A similar formula was proposed by McKee (1989), 

.max 

Further support is provided by rotation. For a core 
with specific angular momentum j, the maximum stable 

1 /2 

mass is given by M^ax ~ [M^ag.max + {4:.8csj/Gf] 
(Tomisaka et al, 1989). The dynamical runaway collapse 
begins when the core mass exceeds this maximum stable 
mass (magnetically supercritical cloud). Quasi-static equi- 
librium configurations exist for cores less massive than the 
maximum stable mass. The evolution of these subcritical 
cores is controlled by the processes of ambipolar diffusion 
and magnetic braking, both of which have longer timescales 
than the gravitational free fall. As the core contracts, the 
density grows and, when n ^ lO^^cm"^, the magnetic field 
is effectively decoupled from the gas. At these densities, 
Joule dissipation becomes important and particle drifts are 
quahtatively different from ambipolar diffusion {Nakano et 
al, 2002). 

The magnetic field is also responsible for the transfer of 
angular momentum in magnetized rotating cores, by a pro- 
cess called magnetic braking. Magnetic breaking is caused 
by the azimuthal component of the Lorentz force {j x B)^. 
In the evolution of subcritical cores, the magnetic braking 
is important during the quasi-static contraction phase con- 
trolled by the ambipolar diffusion (Basu and Mouschovias, 
1994). In the dynamical runaway collapse, the rotational 
speed is smaller than the inflow speed {Tomisaka, 2000) 

3.4. Outflows 

The magnetic field generates an outflow, by which star 
forming gas loses its angular momentum and accretes onto a 
protostar. Magneto-hydrodynamical simulations of the con- 
traction of molecular cores {Tomisaka, 1998, 2000, 2002; 
Allen et al, 2003; Banerjee and Pudritz, 2006) have shown 
that after the formation of the first core, the gas rotates 
around the core, a toroidal magnetic field is induced and 
magnetic torques transfer angular momentum from the disk 
midplane to the surface. Outflows are accelerated in two 
ways: i) The gas which has received enough angular mo- 
mentum compared with the gravity is ejected by the ex- 
cess centrifugal force (magnetocentrifugal wind mecha- 
nism; Blandford and Payne, 1982). ii) In a core with a 
weak magnetic field, the magnetic pressure gradient of the 
toroidal magnetic field accelerates the gas and an outflow is 
formed in the direction perpendicular to the disk. 

Axisymmetric 2D simulations show that i) at least 10% 
of the accreted mass is ejected; ii) the angular momentum is 
reduced to a factor 10~^ of the value of the parent cloud at 
the age of ~ 7000yr from the core formation. This resolves 
the angular momentum problem {Tomisaka, 2000). 7000 yr 
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from the first core formation, the mass of the core reaches 
~ O.IM© and the outflow extends to a distance from the 
core of ~ 2000AU with a speed of ~ 2km (Tomisaka, 
2002). If the accretion continues and the core mass grows 
to one solar mass, the outflow expands and its speed is fur- 
ther increased. It should be noted that the outflow refers to 
the physics of the first core collapse only; the energetics of 
outflows during the second core collapse phase are yet to be 
determined. 

3.5. Radiative Transfer in Multi-Dimensions 

Radiation transport has been shown to play a significant 
role in the outcome of fragmentation into binary and multi- 
ple systems (Boss etal., 2000; Whitehouse and Bate, 2005) 
and in limiting the largest stellar mass in 2D (Yorke and 
Sonnhalter, 2002) and 3D simulations (Edgar and Clarke, 
2004; Krumholz et al, 2005c). The strong dependence of 
the evolution of isothermal and nonisothermal cloud models 
on the handling of the cloud's thermodynamics implies that 
collapse calculations must treat the thermodynamics accu- 
rately in order to obtain the correct solution {Boss et al, 
2000). Because of the great computational burden imposed 
by solving the mean intensity equation in the Eddington ap- 
proximation (the computational time is increased by a fac- 
tor of 10 or more) it is tempting to sidestep the Edding- 
ton approximation solution altogether and employ a sim- 
ple barotropic prescription (e.g.. Boss, 19SI; Bonnell, 1994; 
Bonnell and Bate, 1994a,; Burkert et al., 1997; Klein et al., 
1999). However, Boss et al. (2000) showed that a simple 
barotropic approximation is insufficient and radiative trans- 
fer must be used. We discuss the various methods of radia- 
tion transport in section 3.b.6. 

4. METHODOLOGY OF NUMERICAL SIMULA- 
TIONS 

4.1. Complexity of the Problem of Low Mass Star For- 
mation 

The computational challenge for simulations of low 
mass star formation is that star formation occurs in clouds 
over a huge dynamic range of spatial scales, with different 
physical mechanisms being important on different scales. 
The gas densities in these clouds also varies over many 
orders of magnitude as a result of supersonic turbulence 
and gravitational collapse. Gravity, turbulence, radiation 
and magnetic fields all contribute to the star formation pro- 
cess. Thus the numerical problem is multi-scale, multi- 
physics and highly non-linear. To develop a feel for the 
range of scale a simulation must cover, we can consider 
the internal structure of GMCs as hierarchical, consisting 
of smaller subunits within larger ones (Elmegreen and Fal- 
garone, 1996). GMCs vary in size from 20 to 100 pc, in 
density from 50 to 100 H2 cm~^ and in mass from 10* to 
10*^ Mn. 



Self-gravity and turbulence are equally important in con- 
trolling the structure and evolution of these clouds. Mag- 
netic fields are hkely to play an important role as weU 
{Heilesetal, \992>;McKee etal, 1993). Embedded within 
the GMCs are dense clumps that may form clusters of stars. 
These clumps are few pc. in size, have masses of a few 
thousand M© and mean densities ~ 10^ H2 cm~^. The 
clumps contain dense cores with radii ^ O.lpc, densities 
IC'-IO^ Ho cm and masses ranging from 1 to several 
M0. These cores likely form individual stars or low or- 
der multiple systems. The role of turbulence and magnetic 
fields in the fragmentation of molecular clouds has been in- 
vestigated by 3D numerical simulations (e.g., Padoan and 
Nordlund, 1999; Ostriker et al., 1999; Ballesteros-Paredes, 
2003; Mac Low and Klessen, 2003; Nordlund and Padoan, 
2003). 

A simulation that starts from a region of a turbulent 
molecular clouds (R ^ few pc.) and evolves through the 
isothermal core collapse into the formation of the first hy- 
drostatic core at densities of 10^'^ H2 cm^'^ requires an ac- 
curate calculation across 10 orders of magnitude in density 
and 4-5 orders of magnitude in spatial scale. To resolve 
100 AU separation binaries, one needs a resolution of about 
10 AU. To follow the collapse all the way to an actual star 
would require a further 10 orders of magnitude increase in 
density and 2-3 more orders of magnitude in spatial scale. 
Such extraordinary computational demands rule out fixed 
grid simulations entirely and can be addressed only with 
accurate AMR or SPH approaches. 

4.2. Smooth Particle Hydrodynamics 

The description of the gravitational collapse requires a 
large dynamic range of spatial resolution. An efficient way 
to achieve this is to use Lagrangian methods. Smoothed 
particle hydrodynamics (SPH) is a fully Lagrangian parti- 
cle method designed to describe compressible fluid dynam- 
ics. This method is economical in handling hydrodynamic 
problems that have large, almost empty regions. A vari- 
ety of astrophysical problems including star formation have 
been studied with SPH, because of its simplicity in pro- 
gramming two- and three-dimensional codes and its ver- 
satility to incorporate self-gravity. A broad discussion of 
the method can be found in a review by Monaghan (1994). 
An advantage of SPH is its conservation property; SPH is 
GaUlean invariant and, in contrast to grid-based methods, 
conserves both linear and angular momentum simultane- 
ously. The method to conserve the total energy within a 
computer round-off-error is explained in Inutsuka (2002). 
In order to further increase the dynamic range of spatial 
resolution, Kitsionas and Whitworth (2002) introduced par- 
ticle splitting, which is an adaptive approach in SPH. 

The "standard" SPH formalism adopts artificial viscos- 
ity that mimics the classical von-Neumann Richtmeyer vis- 
cosity. This tends to give poor performance in the descrip- 
tion of strong shocks. In two- or three-dimensional calcu- 
lations of colhding streams, standard SPH particles often 
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penetrate into the opposite side. This unphysical effect can 
be partially eliminated by the so called XSPH prescription 
{Monaghan, 1989), which does not introduce the (required) 
additional dissipation, but results in additional dispersion 
of the waves. As a more efficient method for handling 
strong shocks in the SPH framework, the so called "Go- 
dunov SPH" was proposed by Inutsuka (2002), who imple- 
mented the exact Riemann solver in the strictly conservative 
particle method. This was used in the simulation of the col- 
lapse and fragmentation of self-gravitating cores (Tsuribe 
and Inutsuka, 1999a; Cha and Whitworth, 2003a,b). 

The implementation of self-gravity in SPH is relatively 
easy and one can use various acceleration methods, such as 
Tree-Codes, and special purpose processors (e.g., GRAPE 
board). The flux-limited diffusion radiative transfer was re- 
cently incorporated in SPH by Whitehouse and Bate (2004), 
Whitehouse et al. (2005) and Bastien, Cha, and Viau 
(2004). 

Several groups are now using "sink particles" to follow 
the subsequent evolution even after protostars are formed 
{Bate et al, 1995). This is a prescription to continue the 
calculations without resolving processes of extremely short 
timescale around stellar objects. Krumholz. et al, (2004) 
have introduced sink particles for the first time into Eulerian 
grid-based methods and in particular for AMR. 

4.3. Fixed-Grid Hydrodynamics 

Since the time of its introduction, the numerical code 
of choice for supersonic hydrodynamic turbulence has been 
the Piecewise ParaboUc Method (http://www.lcse.umn.edu/) 
(PPM) of Colella and Woodward (1984). PPM is based on 
a Rieman solver (the discretized approximation to the solu- 
tion is locally advanced analytically) with a third order ac- 
curate reconstruction scheme, which allows an accurate and 
stable treatment of strong shocks, while maintaining nu- 
merical viscosity to a minimum away from discontinuities. 
Because the physical viscosity is not explicitly computed 
(PPM solves the Euler equations), large scale PPM flows 
are characterized by a very large effective Reynolds num- 
ber {Porter and Woodward, 1994). Direct numerical sim- 
ulations (DNS) of the Navier-Stokes equation, where the 
physical viscosity is explicitly computed, require a Unear 
numerical resolution four times larger than PPM to achieve 
the same wave-number extension of the inertial range of 
turbulence as PPM {Sytine et al., 2000). From this point of 
view, therefore, PPM has a significant advantage over DNS 
codes. Furthermore, DNS codes are generally designed 
for incompressible turbulence, and hence of limited use for 
simulations of the ISM. 

Codes based on straightforward staggered mesh finite 
difference methods, rather than Rieman solvers, have also 
been used in applications to star formation and interstellar 
turbulence, such as the Zeus code (http://cosmos.ucsd.edu/) 
{Stone and Norman, 1992a,b) and the Stagger Code 
(www.astro.ku.dk/StaggerCode/) {Nordlund and Gals- 
gaard, 1995; Gudiksen and Nordlund, 2005). Finite dif- 



ference codes address the problem of supersonic turbulence 
with the introduction of localized numerical viscosity to 
stabiUze the shocks while keeping viscosity as low as pos- 
sible away from shocks. The main advantage of this type 
of code, compared with Rieman solvers, is their flexibil- 
ity in incorporating new physics and their computational 
efficiency. 

Fixed-grid codes cannot achieve the dynamical range re- 
quired by problems involving the gravitational collapse of 
protostellar cores. Such problems are better addressed with 
particle methods such as SPH, or by generalizing the meth- 
ods used for fixed-grid codes into AMR schemes. The main 
advantage of large fixed-grid experiments is their abiUty to 
simulate the physics of turbulent flows. As supersonic tur- 
bulence is beheved to play a crucial role in the initial frag- 
mentation of star-forming clouds, fixed-grid codes may still 
be the method of choice to generate realistic large scale ini- 
tial conditions for the collapse of protostellar cores. Recent 
attempts of simulating supersonic turbulence with AMR 
methods are promising {Kritsuk et al., 2006), but may be 
truly advantageous only at a resolution above ~ 1, 000^. 
SPH simulations to date have resolution far too small for 
the task, as commented above, and have not been used so 
far as an alternative method to investigate the physics of 
turbulence. 

4.4. Adaptive Mesh Refinement Hydrodynamics and 
Nested Grids 

The adaptive mesh refinement (AMR) scheme utiUzes 
underlying rectangular grids at different levels of resolu- 
tion. Linear resolution varies by integral refinement factors 
between levels, and a given grid is always fully contained 
within one at the next coarser level (excluding the coars- 
est grid). The origin of the method stems from the seminal 
work of Berger and Oliger (1984) and Berger and Collela 
(1989). The AMR method dynamically resizes and reposi- 
tions these grids and inserts new, finer ones within them ac- 
cording to adjustable refinement criteria, such as the numer- 
ical Jean's condition {Truelove et al., 1997). Fine grids are 
automatically removed as flow conditions require less reso- 
lution. During the course of the calculation, some pointwise 
measure of the error is computed at frequent intervals, typ- 
ically every other time step. At those times, the cells that 
are identified are covered by a relatively small number of 
rectangular patches, which are refined by some even inte- 
ger factor. 

Refinement is in both time and space, so that the calcu- 
lation on the refined grids is computed at the same Courant 
number as that on the coarse grid. The finite difference ap- 
proximations on each level of refinement are in conserva- 
tion form, as is the coupling at the interface between grids 
at different levels of refinement. AMR has three substantial 
advantages over standard SPH. Combined with high order 
Godunov methods, AMR achieves a much higher resolution 
of shocks. This is important in obtaining accuracy in super- 
sonic turbulent flows in star forming clumps and cores and 
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in accretion shocks onto forming protostars. AMR allows 
high resolution at all points in the flow as dictated by the 
physics. Unlike SPH, where particles are taken away from 
low density regions, where accuracy is lost, and concen- 
trated into high density regions, AMR maintains high accu- 
racy everywhere. An important consequence of this is that 
if SPH were to maintain the same comparable resolution 
as AMR everywhere in the flow, it would be prohibitively 
expensive. AMR is based on fixed Eulerian grids and thus 
can take advantage of sophisticated algorithms to incorpo- 
rate magnetic fields and radiative transfer. This is far more 
difficult in a particle-based scheme. AMR was first intro- 
duced into astrophysics by Klein etal. (1990, 1994) and has 
been used extensively both in low mass and high mass star 
formation simulations (Truelove et at., 1998; Klein, 1999; 
Klein etal, 2000, 2003, 2004a; Krumholz etal, 2005c). 

An advantage of SPH over Cartesian grid based AMR is 
that for pure hydrodynamics it can conserve both hnear and 
angular momentum simultaneously to within round-off er- 
rors whereas Cartesian grid based AMR cannot. However, 
if one uses a cylindrical or spherical coordinate system for 
the simulation of protostellar disks for instance, then grid 
based AMR conserves total angular momentum to round- 
off. We point out, however, that these statements apply only 
to pure hydrodynamics. Once forces such as gravity are in- 
cluded, the situation becomes worse and both grid codes 
that solve the Poisson equation or SPH codes that use tree- 
type acceleration or grid based methods for gravity lose the 
conservation property for total linear momentum as well. 

There are several ways to implement AMR. They can 
be broadly divided into two categories: Meshes with fixed 
number of cells, such as in Lagrangian or rezoning ap- 
proaches, and meshes with variable number of cells, such 
as unstructured finite elements, structured cell-by-cell and 
structured sub-grid blocks. For various reasons the most 
widely adopted approaches in astrophysics are structured 
sub-grid blocks and cell-by-cell. The first was developed by 
Berger and Oliger (1984) and Berger and Collela (1989). 
It is used in the AMR code ORION developed by Klein 
and collaborators {Klein, 1999; Crockett et al, 2005) and 
in the community code ENZO {Norman and Bryan, 1999). 
The cell-by-cell approach such as in PARAMESH (Mac- 
Neice et al., 2000) is used in the community code Flash 
{Banerjee et al, 2004). A hybrid approach is used in the 
code NIRVANA {Ziegler, 2005). The cell-by-cell method 
has the advantages of flexible and efficient refinement pat- 
terns and low memory overhead and the disadvantages of 
expensive interpolation and derivation formulas and large 
tree data structures. The sub-grid block method is more ef- 
ficient and more suitable for shock capturing schemes than 
the cell-by-cell method, at the price of some memory over- 
head. 

Finally, nested grids consisting of concentric hierarchi- 
cal rectangular subgrids can also be very effective for prob- 
lems of well defined geometry {Yorke et al, 1993). These 
methods are advantageous for tracing the non-homologous 
runaway collapse of an initially symmetrical cloud in which 



the coordinates of a future dense region are known in ad- 
vance {Tomisaka, 1998). The finest subgrid is added dy- 
namically when spatial resolution is needed as in AMR 
methods. 

4.5 Approaches for Magneto-Hydrodynamics 

Since strong shocks often appear in the astrophysical 
phenomena, a shock-capturing scheme is needed also in 
MHD. Upwind schemes based on the Riemann solver are 
used as the MHD engine. Schemes well known in hy- 
drosimulations, such as Roe's approximate Riemann solver 
{Brio and Wu, 1988; Ryu and Jones, 1995; Nakajima and 
Hanawa, 1996), piecewise parabolic method (PPM; Dai 
and Woodward, 1994), are also applicable to MHD. 

Special attention should be paid to guarantee dxvB = 
in MHD simulations. To ensure that the divergence of 
MaxweU stress tensor = -{l/4TT)BiBj + {l/8TT)B'^Sij 
gives the Lorentz force, the first term of right-hand side 
dj (BiBj) rnxxstheequalto BjdjBi. This requires BidjBj = 
and means that a fictitious force appears along the mag- 
netic field if the condition of divergence-free is broken. 
The divergence of the magnetic field amplifies the insta- 
bility of the solution even for a linear wave. Thus, it is 
necessary for the MHD scheme to keep the divergence of 
the magnetic field zero within a round-off error or at least 
small enough. This divergence-free nature should be satis- 
fied for the boundaries of subgrids in AMR and nested grid 
schemes. 

One solution is based on "constrained transport (CT)" 
{Evans andHawley 1988), in which the staggered colloca- 
tion of the components of magnetic field on the cell faces 
makes the numerical divergence vanish exactly. In the stag- 
gered collocation, the electric field —vxBof the induc- 
tion equation di,B = V x (v x B) is evaluated on the edge 
of the cell-face and the line integral along the edge gives 
the time difference of a component of the magnetic field. 
Note that the electric field on one edge appears twice to 
complete the induction equation. To guarantee a vanishing 
divergence of the magnetic field, CT requires the two eval- 
uations to coincide with each other. 

To utilize the Godunov-type Riemann solver in the con- 
text of CT, Balsara and Spicer (1999) proposed a scheme 
as follows: (1) face-centered magnetic field is interpolated 
to the cell center; (2) from the cell-centered variables, the 
numerical flux at the cell face is obtained using a Riemann 
solver; (3) the flux is interpolated to the edge of the cell-face 
and the electric field in the induction equation is obtained; 
(4) new face-centered magnetic field is obtained from the 
induction equation. Variants of this method are widely used 
[see also Ryu et al, (1998) and Ziegler, (2004)]. 

Avoiding staggered collocation of the magnetic field re- 
quires divergence cleaning. In this case, divergence clean- 
ing is realized by replacing the magnetic field every step 
as S"^"' = B - V$, where V2$ = divB (Hodge pro- 
jection), or by solving a diffusion equation for div B as 
dtB = r/V(V • B). The former is combined with pure 
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Godunov-type Riemann solvers using only cell-centered 
variables (Ryu et al, 1995). Crockett et al. (2005) reported 
that the divergence cleaning of the face-centered magnetic 
field appearing in the numerical flux based on an unsplit, 
cell-centered Godunov scheme improves its accuracy and 
stability. 

Powell et al. (1999) proposed a different formalism, in 
which div B term is kept in the MHD equations as a source 
(e.g., the Lorentz force {W xB)xB/A-k gives an extra term 
related to div B as — SV • B besides the Maxwell stress 
tensor term —Tij, in the equation for momentum density.) 
In this formalism, div B is not amplified but advected along 
the flow. Comparison between various methods is found in 
Toth (2000), Dedner et al. (2002), Balsara and Kim (2004) 
and Crockett et al. (2005). 

There have been attempts to solve the induction equa- 
tion with SPH methods (e.g., Stellingwerf and Peterkin, 
1994; Byleveld and Pong racic, 1996; Price and Monaghan, 
2004a,b,c, 2005). A major obstacle is an instability that de- 
velops when the momentum and energy equations are writ- 
ten in conservation form. As a result, the equations must be 
written in a way that does not conserve momentum {Phillips 
and Monaghan, 1985; Morris, 1996), which is a major con- 
cern for the accurate treatment of shocks. Results of re- 
cent tests of the state-of-the-art SPH MHD code {Price and 
Monaghan, 2004c) appear to be rather poor even for very 
mild shocks, and we conclude that MHD with SPH is not 
yet viable for simulations. 

4.6. Approaches for Radiation IVansport 

Several levels of approximation of the radiation trans- 
port in star formation simulations can be used and details 
of various methods can be found in Mihalas and Mihalas 
(1984) and Castor (2004). Here we briefly describe these 
methods and point out their strengths and weaknesses. Al- 
though modem simulations using radiative transfer are still 
at an early stage, we include methods that hold promise for 
the future that will circumvent the weaknesses of more ap- 
proximate approaches currently in place. 

The simplest improvement beyond a barotropic stiffened 
EOS is the diffusion approximation which pertains to the 
limit in which radiation can be treated as an ideal fluid with 
small corrections. The approximation holds when the pho- 
ton mean free path is small compared with other length 
scales. The combined energy equation for the gas and radi- 
ation results in an implicit non-linear diffusion equation for 
the temperature. The weakness of the diffusion approxima- 
tion is that it is strictly applicable to optically thick regimes 
and performs poorly in optically thin regions. This can be 
severe in optically thin regions of an inhomogeneous turbu- 
lent core or in the optically thin atmosphere surrounding a 
developing protostar. 

The next level of approximation is the Eddington ap- 
proximation (Boss and Myhill, 1992; Boss et al., 2000). 
It can be shown that the diffusion approximation leads di- 
rectly to Eddington's approximation = \EvI, where Pv 



is the pressure tensor moment of the specific intensity of ra- 
diation, Ey is the scalar energy density of radiation and 1^ 
is the isotropic identity tensor. This approximation, coupled 
with dropping the time dependent term in the 2nd moment 
equation of transfer results in a combined parabolic 2nd or- 
der time dependent diffusion equation for the energy den- 
sity of the radiation field. This formulation of the Eddington 
approximation is used in Boss et al. (2000). The approx- 
imation results in a loss of the finite propagation speed of 
light c and a loss of the radiation momentum density, thus 
there is an error in the total momentum budget. In optically 
thin regions, the radiation flux can increase without limit. 
As with all diffusion like methods, this approach also suf- 
fers from shadow effects whereby radiation will tend to fill 
in behind optically thick structures and may lead to unphys- 
ical heating. 

An alternative approach is the Flux Limited Diffusion 
(FLD), that modifies the Eddington approximation, and 
compensates for the errors made in dropping the time de- 
pendent flux term by including a correction factor in the 
diffusion coefficient for the radiation flux. This correction 
factor, called a flux limiter, is in general a tensor and has 
the property that the flux goes to the diffusion limit at large 
optical depth and it correctly limits the flux to be no larger 
than cE in the optically thin regime. This improvement over 
the Eddington approximation has been used by Klein et al. 
(2004c) for the simulation of both low mass and high mass 
star formation. The resulting sparse matrices introduced by 
the diffusion like terms are solved by multi-grid iterative 
methods in an AMR framework. The flux-limiting correc- 
tion can cause errors of order 20in the flux-limiter (or the 
flux), similar to the errors of the Eddington approximation 
of 20% in the Eddington factor at t = in the Milne prob- 
lem (Castor 2004). The FLD metiiod's also suffer from 
shadow effects which can be severe. 

The next level of approximation, the variable Eddington 
tensor method, removes many of the inaccuracies of the Ed- 
dington approximation and the flux limiter modification. It 
was first formulated in multi-dimensions for astrophysical 
problems by Dykema et al. (1996). In essence, if the pre- 
cise ratio of the pressure tensor to the energy density were 
included as an ad hoc multipUer in the Eddington approxi- 
mation equations they would represent an exact closure of 
the system. The tensor ratio is obtained iteratively from ei- 
ther an auxihary solution of the exact transport equation for 
the specific intensity or using an approximate analytic rep- 
resentation of the tensor. This method holds promise for 
future simulations, but has yet to be used in star formation. 

The final two approaches, which are highly accurate and 
deal with the angle dependent transport equation directly, 
are Sn methods and Monte Carlo methods. They have not 
yet been developed for simulations in star formation be- 
cause the cost in 3D is prohibitive. The Sn method is a 
short characteristic method in which a bundle of rays is cre- 
ated at every mesh point and are extended in the upwind di- 
rection only as far as the next spatial cell. The main problem 
is in finding the efficient angle set to represent the radiation 
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field in 2 or 3 dimensions {Castor, 2004). Finally, one might 
consider Monte Carlo methods to solve the transport equa- 
tion. Although simple to implement (its great advantage), 
this method suffers from needing a vast number of oper- 
ations per timestep to get accurate statistics in following 
the particles used to track the radiation field. Both of these 
methods will avoid shadow effects and may be necessary 
to accurately treat optically thick inhomogeneous structures 
that form in accretion flows onto protostars. 

Radiative transfer implementations have recently been 
developed also for SPH methods, based on the flux-limited 
diffusion (Whitehouse et ai, 2005) or the Monte Carlo 
method {Stamatello and Whitworth, 2003, 2005). 

4.7. Comparison of Computational Methods 

Based on the physical processes and numerical method- 
ologies discussed in the previous sections we can compare 
numerical schemes according to their ability to handle the 
following problems both accurately and efficiently: (a) tur- 
bulence, (b) strong shocks, (c) self-gravity, (d) magnetic 
fields, (e) radiative transfer. 

The standard SPH method has been successful with (c) 
and implementations of (e) have been recently developed in 
the flux-limited diffusion approximation and with a Monte 
Carlo method. It does not include (d), it is well known to be 
inadequate for (b) and has had virtually no apphcations to 
(a) to date. As any Lagrangian particle methods, SPH pro- 
vides good resolution in high density regions, but very poor 
in low density ones. The Godunov SPH method improves 
the standard SPH codes because of its ability to address (b), 
but does not provide a solution to (d) and is untested for 
(a) as well. Although MHD is currently under development 
in SPH, results of standard MHD tests with a state-of-the- 
art code show the need for significant improvements even 
in the case of very mild shocks. Current applications of 
SPH should therefore be limited to non-MHD problems and 
the accuracy and performance of SPH with turbulent flows 
must be thoroughly tested. 

In hydrodynamical problems, grid-based methods such 
as MUSCL (van Leer, 1979) and PPM (Colella and Wood- 
ward, 1984) use exact Riemann Solvers to construct the nu- 
merical fluxes and provide very accurate description in as- 
trophysical flows with strong shocks (b). They have also 
been thoroughly tested with compressible turbulent flows, 
and MHD versions have been developed that can address 
both (d) and (e). Traditional finite-difference grid-based 
schemes may still be useful, because the best of them can 
also accurately address (a), (b), (d) and (e), at a lower cost 
of code development and computer resources. Point (c) can 
also be efficiently dealt with by grid-based codes thanks 
to AMR methods. However, the development of AMR 
schemes that satisfy (c) and at the same time (d) has be- 
gun only recently. These schemes exist and have been suc- 
cessfully tested, but it is unclear at present which approach 
will provide the best trade off between accuracy and perfor- 
mance. 



The constrained transport method appears to be the ideal 
one to guarantee the V ■ -B = condition. Various schemes 
have been proposed even in the category of Godunov-type 
methods with a linearized Riemann Solver. An exact MHD 
Riemann Solver would be more adequate to handle strong 
shocks, but it is not available yet. In MHD we have to 
solve seven characteristics even in one-dimensional prob- 
lems, which hinders an efficient construction of numerical 
fluxes based on the non-linear waves. Furthermore, the dis- 
covery of the existence of the MHD intermediate shocks 
{Brio and Wu, 1988) brings an additional difficulty in the 
categorization and prediction of the emerging non-Unear 
waves. Among possible solutions, a linearized Riemann 
Solver with artificial viscosity may still be a useful option. 

The Godunov MHD code of Crockett et al, (2005) has 
been merged with the AMR RHD code of Klein et al. 
(2004a,b) into the first fully developed AMR magneto- 
radiation-hydrodynamic code (ORION) to be used in simu- 
lations of star formation capable of addressing (a) through 
(e). 

5. RECENT SIMULATIONS AND CONFRONTA- 
TION WITH THE OBSERVATIONS 

5.1. Turbulent Fragmentation of Molecular Clouds 

The fragmentation of molecular clouds is the result of 
a complex interaction of supersonic turbulence with grav- 
ity and magnetic fields. Supersonic turbulent flows gener- 
ate nonUnear density enhancements through a complex net- 
work of interacting shocks. Some density enhancements 
are massive and dense enough to undergo gravitational in- 
stability and collapse into stars. A fundamental problem 
with our understanding of star formation is that the physics 
of turbulence is not fully understood. In order to investi- 
gate the process of turbulent fragmentation we must rely on 
large numerical simulations that barely resolve the scale- 
free nature of interstellar turbulent flows. If the scaling 
laws of turbulence play a role in the star formation process, 
we cannot accurately test their effects numerically, unless 
those scaling laws are well reproduced in the simulations. 
Ballesteros-Paredes et al. (2006) have tested the idea that 
the power spectrum of turbulence determines the Salpeter 
IMF {Padoan and Nordlund, 2002). However, their simula- 
tions do not generate a turbulence inertial range, due to the 
combined effect of low resolution and large numerical dif- 
fusivity. Their velocity power spectra do not show any ex- 
tended power laws, and they even differ between their grid 
and SPH simulations. As a consequence, such numerical 
simulations fail to reproduce a scale free mass distribution 
of unstable cores. 

Other recent SPH simulations have been used to com- 
pute the stellar mass distribution (e.g., Bonnell et al, 2003; 
Klessen, 2001) and to test the effect of turbulence on star 
formation {Delgado-Donate et al., 2004). Although such 
SPH simulations are ideally suited to follow the collapse of 
individual objects due to their Lagrangian nature, their size 
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is far too small to generate an inertial range of turbulence, 
as discussed in section 3.a.l. 

Despite their limitations, numerical simulations carried 
out over the last few years have given us a good statis- 
tical picture of the process of fragmentation of magne- 
tized supersonic flows. The foUowing are the most im- 
portant results: i) The dissipation time of turbulence is al- 
most independent of the magnetic field strength. The tur- 
bulence decays in approximately one dynamical time in 
both equipartition and super-Alfvenic flows {Padoan and 
Nordlund, 1997; MacLow et al, 1998; Stone et ah, 1998; 
Padoan and Nordlund, 1999; Biskamp and Muller, 2000). 
ii) The velocity power spectrum and structure functions are 
power laws over an inertial range of scales (Boldyrev, 2002; 
Boldyrev et al. 2002a,b; Padoan et al., 2004b). In the 
limit of very large rms Mach number, the turbulent velocity 
power spectrum scales approximately as u'^ (x kr^-^ and 
the velocity structure functions follow the relative scaling 
predicted by Boldyrev (2002). For intermediate levels of 
compressibility, the scaling exponents depend on the rms 
Mach number {Padoan et al, 2004b). iii) The power spec- 
trum of the gas density is a power law over an inertial range 
of scales. Its slope is a function of the rms Mach number of 
the flow and of the average magnetic field strength {Padoan 
et al, 2004b; Beresnyak et al, 2005). iv) With an isother- 
mal equation of state, the probability density function of gas 
density is well approximated by a Log-Normal distribution 
{Vazquez-Semadeni, 1994; Padoan etal, 1997; Scalo etal, 
1998; Passot and Vazquez-Semadeni, 1998; Nordlund and 
Padoan, 1999; O striker et ah, 1999; Wada and Norman, 
2001), with the dispersion of linear density proportional to 
the rms Mach number of the flow {Padoan et al, 1997; 
Nordlund and Padoan, 1999; O striker et ah, 1999; Li et 
al, 2004). v) If the kinetic energy exceeds the magnetic 
energy, the distribution of the magnetic field strength, B, is 
very intermittent and is correlated with the gas density, n. 
The scatter plot of B versus n shows a very large dispersion 
and a well defined power law upper envelope {Padoan and 
Nordlund, 1999; Ostriker et al, 2001). If the kinetic en- 
ergy is comparable to the magnetic energy, strong density 
enhancements are still possible in the direction of the mag- 
netic field, but fluctuations of B are small and independent 
of n. vi) The flow velocity is correlated to the gas density. 
Because density is increased by shocks, where the velocity 
is dissipated, dense filaments and cores have lower velocity 
than the low density gas {Padoan et al, 2001b). vii) The 
mass distribution of gravitationally unstable turbulent den- 
sity peaks is very close to the stellar IMF and follows the 
analytical model of Padoan and Nordlund (2002) {Li et al , 
2004; Tilley and Pudritz, 2004). 

There is now substantial observational evidence indicat- 
ing that these main properties of supersonic MHD turbu- 
lence are indeed found in molecular clouds. The compar- 
ison of numerical simulations of turbulence with observa- 
tional data was pioneered by Falgarone et al (1994) and 
continued by many others (e.g., Padoan et al. 1998, 1999; 
Padoan et al. 2001a,b; Ostriker et al, 2001; Ballesteros- 



Paredes and Mac Low, 2002; Ossenkopf, 2002; Padoan et 
al, 2004a; Gammie et al, 2003; Klessen et al., 2005; Es- 
quivel andLazarian, 2005). 

However, with the exception of several works by Padoan 
et al. and the work of Ossenkopf (2002), where post- 
processed three dimensional non-LTE radiative transfer cal- 
culations were carried out, all other studies are based on a 
simple comparison of density and velocity fields in the sim- 
ulations with the observed quantities. Some recent works 
addressing the comparison of turbulence simulations and 
observational data include studies of velocity scaling, show- 
ing that molecular cloud turbulence is driven on large scale 
(e.g., Ossenkopf and Mac Low, 2002; Heyer and Brunt, 
2004) and studies of core properties, showing that turbu- 
lent flows naturally generate dense cores with shapes, inter- 
nal turbulence, rotation velocity and magnetic field strength 
consistent with the observations (e.g., Padoan and Nord- 
lund, 1999; Gammie et al, 2003; Tilley and Pudritz, 2004; 
Tilley and Pudritz, in preparation; Li et al, 2004). 

5.2. Collapse and Fragmentation of Molecular Cloud 
Cores into Low-Mass Stars 

Over the past several years two dominant models of how 
stars acquire their final mass have emerged: Direct Gravi- 
tational Collapse and Competitive Accretion. In both the- 
ories, a star initially forms when gravitational bound gas 
collapses. In the gravitational collapse scenario, after a 
protostar has consumed or expelled all the gas in its initial 
core, it may continue accreting from the parent clump, but it 
will not significantly alter its mass {McKee and Tan, 2003; 
Padoan et al, 2005; Krumholz et al, 2005b). Competitive 
accretion, in contrast, requires that the amount accreted af- 
ter consuming the initial core be substantially larger than 
the protostellar mass. 

Krumholz et al (2005a) define fm = int:tdyn/ "i* as the 
fractional change in mass that a protostar of mass m* un- 
dergoes each dynamical time tdyn of its parent clump, start- 
ing after the initial core has been consumed by the accret- 
ing protostar. Gravitational collapse theory suggests that 
Im 1, while competitive accretion requires f„i ^ 1. In 
recent work examining the plausibility of competitive ac- 
cretion, Krumholz et al (2005a) considered two possible 
scenarios. In the first scenario, the gas the protostar is ac- 
creting is not accumulated into bound structures on scales 
smaller than the entire clump. For unbound gas, self-gravity 
may be neglected and the entire problem can be treated as 
Bondi-Hoyle accretion in a turbulent medium of non-self- 
gravitating gas onto a point particle. In a companion paper 
{Krumholz et al, 2006) they develop the theory for Bondi- 
Hoyle accretion in a turbulent medium. Using this theory 
they derive the accretion rate for such a turbulent medium 
and they confirm their theory with detailed, high resolution, 
converged AMR simulations. By using this accretion rate 
and the definition of the virial parameter, avir = vir/^, 
and the dynamical time, t^yn = R/<y, where a is the veloc- 
ity dispersion in the gas, they show that the accretion of un- 
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bound gas gives /m-BH = (14.4, 3.08^)(/)BHavir ^(^) 
for a (spherical, filamentary) star-forming region, where 
0BH represents the effects of turbulence {Krumholz et al, 
2005a). 

From this it is clear that competitive accretion is most ef- 
fective in low mass clumps with virial parameters ctvir ^ 1- 
They then examined the observed properties of a large range 
of star forming regions spanning both low mass and high 
mass stars and computed the properties for each region 
yielding c^vir, f^BH and /m-BH- In virtually every region 
examined, the virial parameter ctvir ^ 1 and /m-BH •C 1. 
Thus none of the star-forming regions are consistent with 
competitive accretion, but all are consistent with direct 
gravitational collapse. Edgar and Clarke, (2004) exam- 
ined Bondi-Hoyle accretion onto stars including radiation 
pressure effects and found that radiation pressure halts fur- 
ther accretion around stars more massive than ^ 1GM0. It 
is important to point out that while this result may be true 
for accretion of unbound gas onto a point particle, it is not 
true for global collapse and accretion from a bound core 
as shown in more reahstic full 3D radiation hydrodynamics 
simulations by Krumholz et al. (2005c), who form massive 
stars with M ^ 4OAf0. If Edgar and Clarke's results do 
hold for Bondi-Hoyle accretion (but not accretion from a 
core), then the only way for massive stars to grow in com- 
petitive accretion is by direct collisions. This requires den- 
sities of 10^ pc""^, 3 orders of magnitude larger than any 
observed in the galactic plane. Furthermore, no competitive 
accretion model to date has included the effects of radia- 
tion pressure; a glaring omission if the model is attempting 
to explain high mass stars. It follows that competitive ac- 
cretion is not a viable mechanism for producing the stellar 
IMF. 

In a second possible competitive accretion scenario, 
Krumholz et al. (2005a, supplemental section) examined 
another way that a star could increase its mass by captur- 
ing and accreting other gravitationally bound cores. Their 
theory results in the calculation of /m-cap, the fractional 
change in mass that a protostar undergoes by capturing 
bound cores. As found with /m-BH, all the values are es- 
timated to be three more orders of magnitude below unity 
and again, competitive accretion is found not to work. 

If competitive accretion is clearly not supported by ob- 
servations in any known star forming region, why do the 
simulations (Bonnell et a/., 1998, 2001; Bate et al, 2005) 
almost invariably find competitive accretion to work? Is 
there a fundamental flaw in the methodology used in com- 
petitive accretion scenarios (SPH) or is the problem one of 
physics and initial conditions? As Krumholz et al. (2005a) 
point out, all competitive accretion have virial parameters 
avir 1. Some of the simulations start with ayir ~ 0.01 
as a typical choice (Bonnell et al., 2001a,b; Klessen and 
Burkert, 2000, 2001). For other simulations the virial pa- 
rameter is initially of order unity but decreases to <C 1 in 
a crossing time as turbulence decays (Bonnell et al., 2004; 
Bate et al, 2002a,b, 2003). It is also noteworthy that many 
of the simulations begin with clumps of mass consider- 



ably smaller (M < IOOMq) than that typically one found 
in star forming regions ^ 5OOOM0 (Plume et al, 1997). 
Krumholz et al. (2005a) show that for competitive accre- 
tion to work, avir^M < 10 Mq, but for typical star forming 
regions avir ~ 1 and M « 10^-10^ and the inequahty 
is almost never satisfied. 

One reason why simulations evolve to avir ^ 1 is 
that they omit feedback from star formation. Observa- 
tions by Quillen et al. (2005) show that outflows inject 
enough energy on the scale of a clump to sustain turbu- 
lence thereby keeping the virial parameter from declining 
to values much less than unity. Another possible reason is 
that the simulations consider isolated clumps with too lit- 
tle material. Real clumps are embedded in larger molecu- 
lar clouds where larger scale turbulent motions can cascade 
down to the clump scale preventing the rapid decay of the 
turbulence. 

5.3. 3D Collapse with Radiation 

Radiation transfer is an important element in both low 
mass and high mass star formation. Boss et al. (2000) car- 
ried out the first 3D simulations including radiation trans- 
fer to study the effect of radiation on the formation of 
stars in collapsing molecular cloud cores. Starting from 
cores with Gaussian initial density profiles, this work com- 
pared collapse calculations based on the isothermal and 
barotropic approximations with the more realistic case of 
detailed radiative transfer in the Eddington approximation. 
The use of the isothermal equation of state resulted in a 
collapse leading to a thin isothermal filament (Truelove 
et al, 1997). In the more realistic case with nonisother- 
mal heating using radiative transfer in the Eddington ap- 
proximation, they showed that thermal retardation of the 
collapse caused the formation of a binary protostar sys- 
tem at the same maximum density where the isothermal 
collapse yielded a thin filament. Eventually, the binary 
clumps evolved into a central protostar surrounded by spi- 
ral arms containing two more clumps. The corresponding 
collapse using the barotropic approximation allowed a tran- 
sition from an isothermal optically thin to an optically thick 
flow. It resulted in a transient binary merging into a cen- 
tral object surrounded by spiral arms with no evidence of 
further fragmentation. The barotropic result differs signifi- 
cantly from the Eddington result at the same maximum den- 
sity indicating the importance of detailed radiative transfer 
effects. 

Boss et al. (2000) examined the differences in the use 
of a barotropic stiffened EOS approximation and radiative 
transfer. In the former, the thermal properties of the gas 
are specified solely as a function of density. This implies 
that a single global value of a critical density represents the 
entire calculation and its value depends weakly on the as- 
sumed geometry of the cloud (Inutsuka andMiyama, 1997). 
In 3-D the effective value of this critical density depends 
on the local geometry surrounding a fluid element. In ad- 
dition, whereas the specific entropy of a gas parcel gener- 
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Fig. 1. — Binary outflows. Magnetic field lines, velocity 
field vectors and density distribution on the z = plane. 
Taken from model DL of Machida et al (2005b). 

ally depends on the thermal history of the parcel, the spe- 
cific entropy of the gas using the barotropic approximation 
depends solely on the density. Thus the derived pressures 
used in the momentum equations will differ between a cal- 
culation using a stiffened EOS approximation and a fully 
consistent calculation using radiative transfer As a result, 
the temperature is determined not simply by adiabatic com- 
pression, but by compressional heating in a 3D volume with 
highly variable optical depth. Thus the dependence of the 
temperature on the density cannot be represented with a 
simple barotropic approximation with any great accuracy. 
This causes concern for the validity of current simulations 
of multiple star formation and cluster formation, since es- 
sentially all use either the isothermal or the barotropic ap- 
proximation. Recent work by Whitehouse and Bate (2005) 
examined core collapse with radiation and the adequacy of 
the barotropic approximation as well. 

5.4. The Debate Over Disk Fragmentation 

As pointed out in section 2, most simulations to date, 
performed with SPH, find circumstellar disk fragmenta- 
tion and rapid ejection of BDs in most cases {Bate et al., 
2QQ3;Delgato-Donate et ai, 2004; Goodwin etai, 2004a,b) 
even possibly to an excess of BD and low-mass companions 
(Goodwin et ai, 2004b). Furthermore, most simulations 
are terminated at an arbitrary time, when much of the gas 
is still present; hence the simulations may provide only a 
lower bound to the number of companions that would be 
produced. As discussed in section 2, these simulations are 
increasingly contradicted by recent observations. Goodwin 
and Kroupa (2005) have recently pointed out that obser- 
vations suggest that individual cores produce at most 2-3 



protostars. A possible reason for this problem is the ab- 
sence of a magnetic field and inaccurate thermodynamics 
(i.e. barotropic equation rather than radiative transfer) in 
the SPH simulations. 

However, there is an alternative explanation: The SPH 
simulations are likely not converged. Indeed, recent high 
resolution AMR simulations of the collapse and fragmen- 
tation of turbulent cores {Klein et al., 2004a,b) show that 
a magnetic field and radiative transfer is not required to 
explain the observational results of Goodwin and Kroupa 
(2005). Klein et al, (2004a,b) find that over a broad range 
of turbulent Mach numbers (M ^ 1 — 3) and rotational 
energies {j3 ~ lO""* — 10~^) low order single or binary 
stars are formed through fragmentation of the core, not the 
ensuing disk. Is it then a matter of faith in one numerical 
method or the other? SPH versus AMR? Not really, it is 
rather a matter of testing the convergence of the numerical 
simulations, irrespective of the method. 

Fisher et al. (2006, in preparation) have performed high 
resolution, full convergence studies using the SPH code 
Gadget with the same model and initial conditions as Good- 
win et al. (2004a), except for the initial turbulent seed. 
They confirmed the results of Goodwin et al., (2004a) at the 
same low resolution (only one smoothing kernel per mini- 
mum Jeans mass). But they have also found no convergence 
of either the multiple number of companion protostars pro- 
duced or the time for multiple fragmentation to occur in 
the disk, with up to 32 times the resolution of Goodwin et 
al., (2004a). This suggests that the current SPH simula- 
tions showing high order multiple disk fragmentation may 
be grossly under resolved. If so, the disagreement with the 
observations may not be due only to the absence of the mag- 
netic field and radiation, but also to unconverged numerics. 
Recent grid-based simulations attempting to study fragmen- 
tation in isolated disks (cf. Duriesen this book) indicate that 
for 2D axisymmetric disk studies very high resolution, 256 
cells in the radial direction alone {Ostriker private commu- 
nication), is required to demonstrate that disks do not suffer 
from numerical instability. This level of resolution is far be- 
yond any SPH simulation of cores or clusters to date show- 
ing disk fragmentation, and would be difficult to achieve 
also with AMR simulations starting from globally collaps- 
ing cores. 

At the time of this conference, the debate between SPH 
and AMR with respect to the issue of disk fragmentation 
continues, but it is our strong opinion that all simulations 
using SPH or AMR must demonstrate adequate conver- 
gence to be credible. This should apply also to physical 
systems that may display chaotic behavior, such as the gas 
dynamics in a molecular cloud core. A high sensitivity to 
initial conditions may result in statistical distributions of 
the measured quantities (e.g., number of collapsing objects 
and their formation time) around some mean value. In this 
case, an expensive approach to test numerical convergence 
would consist of running a large number of experiments at 
each resolution and test for the convergence of the statistical 
distributions of the quantities of interest. A less expensive 
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Fig. 2. — Bird's-eye view of the magnetic field lines in the 
fast jet emanating from the central region around a proto- 
star. Only the 17th grid in the nested-grid resistive MHD 
calculations by Machida et al. (2006, in preparation) is 
shown. The density contours and velocity vectors are pro- 
jected on the walls on both sides. 

approach is to measure quantities with a weak sensitivity to 
initial conditions because they already represent the average 
of some statistical distribution, or just because they are even 
more sensitive to the numerical resolution than to the initial 
conditions. For example. Fisher et al. (2006) find that as 
they increase the resolution of their SPH simulations, the 
time of disk fragmentation (when the first object is formed) 
increases monotonically with resolution, a sign of lack of 
convergence in the SPH simulations, rather than a signature 
of chaotic behavior. 

5.5. Star formation in a cloud with Magnetic Field and 
Rotation 

5.5.7. Fragmentation of the First Core. Anon-rotating 
cloud core without magnetic field contracts in a self-similar 
manner to form a first core that is composed of molecu- 
lar hydrogen {Larson, 1969; Masunaga et al, 1998). Re- 
cently, nested grid MHD simulations (Machida et al., 2004, 
2005a) have revealed that i) a rotating magnetized core 
evolves maintaining a ratio of angular speed to magnetic 
field strength at the center Vlc/Bc — const; and ii) J7c and 
Be are well correlated at the first core phase and satisfy 
the "magnetic flux - spin relation" as Oc/0.2^47rGpc + 
i?^/0.36^87rCgpc — 1, using a central density pc and 
isothermal sound speed c^. This is regarded as an appear- 



ance of self- similarity in magnetized rotating cores. 

The fragmentation of the first core is regarded as one 
of the mechanisms to explain close binary systems {Bon- 
nell and Bate, 1994b; Bate, 1998). The magnetic field 
affects the rotational motion (magnetic braking) and thus 
the fragmentation. Whether the magnetic field stabiUzes 
the first core against the fragmentation or not {Machida 
et al, 2005b; Ziegler, 2005) is attracting attention in re- 
lation to the binary formation. As well as non-magnetic 
cores (see §3.a.5), a magnetized first core can fragment 
if it is rotating sufficiently fast, fic ^ G.2(47rG'/9c)^/^ 
{Machida et al., 2005b) and has only a weak magnetic 
field. Be ^ 0.3(87rc^pc)^/^. Simulations show that in- 
creasing the magnetic field strength, fragmentation is sta- 
bihzed by the suppression of rotational motion by mag- 
netic braking (see also Hosting and Whitworth, 2004). 
This is not found by Boss (2002), but his model equation 
is not fully consistent with MHD and does not account 
for magnetic braking. In order to achieve enough rota- 
tion to cause fragmentation, the initial fl-to-B ratio must 
satisfy the condition {Q/B)init > 0.39G^/^Cs ~ 1.7 x 
10-^(c,/0.19kms-i)-VG-^ {Machida et al, 2005b). 

When B and angular momentum J are not parallel to 
each other, the magnetic braking works more efficiently for 
the component of J perpendicular to B. A magnetically 
dominated cloud core whose f2-to--B ratio is less than the 
above critical value forms a disk perpendicular to the mag- 
netic field {Matsumoto and Tomisaka, 2004) and an outflow 
is ejected in the direction of the local magnetic field, even 
if B is not parallel to J. The difference between the local 
field in the vicinity of a protostar and on the cloud scale is 
restricted to < 30 deg for this case. 

If the first core is fragmented into binary or multiple 
cores, each fragment spins and multiple outflows (Fig.l) 
are ejected {Machida et al, 2004, 2005b; Ziegler, 2005; 
Banerjee and Pudritz, 2006), which explains binary out- 
flows {Liseau et al, 2005). 

5.5.2. Second Collapse with Magnetic Field. Once the dis- 
sociation of molecular hydrogen starts at the central region 
of the first core, the second collapse begins. Further cal- 
culation of the evolution to form the second core (i.e., a 
protostar) requires the inclusion of resistivity in the MHD 
description, as the high density reduces the degree of ion- 
ization and the conductivity of the medium. Machida et al 
(2006, in preparation) have adopted parametrized resistivity 
as a function of density in their nested-grid resistive MHD 
code, and extended the calculations of the collapse. Dur- 
ing the isothermal phase, the magnetic Reynolds number is 
a decreasing function of density. If the magnetic Reynolds 
number decreases below unity, the magnetic field is effec- 
tively decoupled from the collapsing gas. However, the 
temperature of the gas becomes sufficiently high (~ lO^K) 
that the magnetic field is re-coupled again with the collaps- 
ing gas. This relatively sudden grab of the field lines tends 
to make a very coUimated fast outflow around the second 
core. Fig. 2 shows a snapshot of the propagation of a fast 
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jet from the protostar. The density distribution on the cross 
section is also shown on the left and right walls. Note that 
bow shocks are clearly seen in the density plot. The re- 
sults of this calculation indicate that a realistic modeling of 
the evolution of temperature and resistivity as a function of 
density is required for a precise description of the jet. 

6.SUMMARY AND FUTURE DIRECTIONS 

6.1. Summary 

Observations of molecular clouds and cores provide a 
wealth of data that are important constraints for initial 
conditions of numerical simulations. Large scale simula- 
tions should be consistent with the turbulent nature of the 
ISM and the observed properties of prestellar cores should 
emerge self-consistently from the simulations. All calcu- 
lations must adhere to the Jeans condition for grid-based 
schemes or a well established equivalent for particle-based 
schemes. It is important to stress that the Jeans condition is 
a necessary but not sufficient condition to guarantee avoid- 
ance of artificial fragmentation. Simulations must be tested 
at more resolved Jeans numbers to establish convergence. 

The importance of turbulence in star formation is now 
well accepted. A very large spatial resolution is required 
to simulate the turbulent fragmentation, barely achieved by 
the largest grid-based simulations. Present SPH simulations 
fall several orders of magnitude below the required spatial 
resolution. It is possible that almost no available simula- 
tions have yet accurately tested the effect of turbulent frag- 
mentation. 

Magnetic fields play a crucial role in the star formation 
process. At this time there are no 3D, self-gravitational 
MHD simulations that have evolved stars from turbulent 
clouds. Grid based schemes such as AMR have devel- 
oped high order accurate PPM and Godunov based MHD 
that can provide accurate solutions across several orders 
of magnitude of collapse. SPH has developed cruder ap- 
proaches to MHD that appear to be rather poor even for very 
mild shocks, but this will hopefully improve. Godunov ap- 
proaches to MHD for SPH may provide increased accuracy. 

Radiation transport has been shown to play a significant 
role in the outcome of fragmentation to binary and multi- 
ple systems. It has been shown that side stepping the issue 
of radiation transport with a barotropic approximation can 
lead to incorrect results. Current AMR calculations have 
implemented radiation transfer in the flux-limited approxi- 
mation and have used this for both low mass and high mass 
star formation. Radiation transfer schemes have recently 
been developed for SPH as well. 

Current star formation simulations are not yet adequate 
to accurately span the many orders of magnitude of den- 
sity and spatial range necessary to account for stars from 
initially turbulent clouds and encompassing all the relevant 
physics. In our opinion, AMR approaches with recent de- 
velopments coupUng radiation transfer and MHD hold out 
the best promise for achieving that goal. SPH, while mak- 
ing significant strides in the last few years, is still faced 



with difficult challenges of accurate handling of turbulence, 
strong shocks, MHD, and radiation transfer. However, we 
anticipate further progress with SPH. 

Current calculations are still not adequate to explain the 
stellar IMF. Of the two dominant theories for the origin 
of stellar masses described as direct gravitational collapse 
and competitive accretion, observations appear to favor the 
first. Recent theoretical work by Knimholz. et al. (2005a) 
has now demonstrated that seed protostars cannot gain mass 
efficiently by competitive accretion processes in observed 
star-forming regions that are approximately in virial bal- 
ance. There is no observational evidence for the existence 
of regions that are far from virial balance, as required by 
competitive accretion models. This suggests that competi- 
tive accretion is not a viable mechanism for producing the 
stellar IMF and that current simulations resulting in com- 
petitive accretion must have initial or boundary conditions 
inconsistent with the observations, or rather neglect some 
crucial physics. Theoretical efforts directed toward the pic- 
ture of direct gravitational collapse set up by turbulent frag- 
mentation appear to be more promising. 

6.2. Future Directions for Numerical Simulation 

The lack of demonstrated convergence for most simula- 
tions in the field presents us with uncalibrated results. We 
strongly emphasize that future simulations should demon- 
strate numerical convergence before detailed comparison 
with observations can be credible. Otherwise there is no 
way to normahze the accuracy of large scale simulations 
and the results of such simulations will not advance our un- 
derstanding of low mass star formation. Convergence tests 
should be always carried out irrespective of the numerical 
approach. A better understanding of the numerical treat- 
ment of disk fragmentation must occur to clear up the cur- 
rent discrepancies between AMR and SPH. 

Future simulations of low mass star formation must en- 
deavor to include MHD and radiation transfer. With the 
development of accurate approaches to these processes, we 
can expect to see simulations become more relevant to ad- 
dressing the observations. For simulations to make a real 
connection to the observations, detailed fine profiles and 
continuum sub-mm and mm maps should be calculated with 
3D radiative transfer codes. Approaches that go beyond the 
Eddington approximation and flux-limited diffusion such as 
Sn transport and Monte Carlo need further development to 
work efficiently in full 3D simulations. They will become 
increasingly important in treating the flow of radiation in 
highly inhomogeneous regions. As future simulations en- 
compass multi-coupled physics, significant progress must 
be made in algorithms that improve the parallel scalabil- 
ity. That is necessary in order to simulate the fuU dynamic 
range of collapse and fragmentation from clouds to stars, 
while capturing all the relevant physics. 
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